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ABSTRACT 

The monitoring of quiescent emission from neutron star transients with accretion outbursts long enough to 
significantly heat the neutron star crust has opened a new vista onto the physics of dense matter. In this paper we 
construct models of the thermal relaxation of the neutron star crust following the end of a protracted accretion 
outburst. We confirm the finding of Shternin et al., that the thermal conductivity of the neutron star crust is 
high, consistent with a low impurity parameter. We describe the basic physics that sets the broken power-law 
form of the cooling lightcurve. The initial power law decay gives a direct measure of the temperature profile, 
and hence the thermal flux during outburst, in the outer crust. The time of the break, at hundreds of days 
post-outburst, corresponds to the thermal time where the solid transitions from a classical to quantum crystal, 
close to neutron drip. We calculate in detail the constraints on the crust parameters of both KS 1731-260 and 
MXB 1659-29 from fitting their cooling lightcurves. Our fits to the lightcurves require that the neutrons do not 
contribute significantly to the heat capacity in the inner crust, and provide evidence in favor of the existence of 
a neutron superfluid throughout the inner crust. Our fits to both sources indicate an impurity parameter of order 
unity in the inner crust. 

Subject headings: dense matter — stars: neutron — X-rays: binaries — X-rays: individual(KS 1731-260, 
MXB 1659-29) 



1. INTRODUCTION 

In 2001 the low mass X-ray binary KS 1731-260 went 
into quiescence ( Wijnands et al. 200 \ \ after a ccreting steadily 
since its first detection about 12 years prior (Su nyaev|1989| l. 
Rutledge et al.| ( |2002[ ) suggested that the neutron star crust 
would be heated out of thermal equilibrium with the core 
during this long outburst, and that monitoring observations 
could detect the thermal relaxation of the crust following 
the cessation of accretion. A regular monitoring program 
of KS 1731-260 with Chandra and XMM ([Wijnands et al.| 
200Tj |Wijnands et al.||2002j |Cackett et al."1|2006| ) has indeed 
detected a steadily decaying luminosity. Another source, 
MXB 1659-29, has also been observed to cool following a 
2.5 yr -long outburst (Wijnands et al. 2003; Wijnands et al.| 
2004), and recent observations show that the cooling ap- 
pears to have halted ( |Cackett e t al. 2008). Recently a third 
source, AX J 1754.2-2754, was observed to turn off ( |Bassa| 



significantly different properties than that of a young isolated 
neutron star, as the accretion lifetime of a LMXB is long 
enough for accreted matter to replace the entire crust. The ac- 
creted crust is heated during outburst by electron captures and 
pycno nuclear reactions, mostly at densitie s close to neutron 
drip ( |Haensel & Zdunik|1990||2lX)3l|2008] ). The temperature 
profile in the crust is determined by how the heat released in 
this "deep cmstal heating" is transported outwards to the sur- 
face or inwards to the neutron star core, where it is radiated as 
neutrinos. 

A particularly uncertain property of accreted crusts is the 
thermal conductivity. The matter entering the top of the crust 
is expected to consist of a mixture of elements produced by 
hydrogen and helium burning at low densities. Schatz et al. 
(1999) found that rp-process burning produced ashes with 
Gimp ~ 100, where the impurity parameter 



et al. 2008) after being in outburst since being detected in 
1999 ( Sakano et al.|20 02 ). This source is likely an ultracom- 
pact binary, and is of particular interest as it may be simi- 
lar to 1H 1905+000, which has an extremely low quiescent 
flux, perhaps indicating a very low neutron star core tem- 
perature (Jonk eTet al.||2006| . Finally, the rapidly accreting 
source XTE J 170 1-462 we nt into quiescence in 2007 after 
1 .6 yr of active accretion ( Homa n et al.|2007| Altamirano 



Gimp = »r o n Yj " i(Zi ~ {Z))1 



(1) 



measures the distribution of nuclide charge numbers, an im- 
portant parameter f or setting the co nductivity of the crust (Itoh 
& Kohyama 1993[ ). Brown (2000 1 suggested that at such large 
values of Gimp, the crust would have an amorphous structure, 
with a low thermal conductivity comparable to the conduc- 
tivity of the liquid state. Recent molecular dynamics calcula- 
et al.|2007] ), and EXO 0748-67 6 entered quiescence i n 2008 tions of the s olidification of such a mixture ( |Horowitz et al.| 



a fter being in outburst fo r 24 yr ( Degena ar et al.|2009] ) 

|Rutledge et al. ( 2002| ) emphasized that observations of the 
crust cooling would offer a new probe of the crust thermal 
properties. The crust of a neutron star in a LMXB may have 
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2007 2008a[ ) show, however, that a regular crystal structure 
does form, but that phase separation between the liquid and 
solid phases or the formation of multiple solid phases can oc- 
cur, changing the effective value of Q lmv in the crust com- 
pared to the rp-process mixture. In addition, nuclear reactions 
in the crust may also act to reduce the effective value of Q lmp 
(Horowitz et al. 2008a| ). Measurement of Gimp in the crust 
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would therefore give an important test of the composition of 
rp-process ashes and their subsequent evolution as they are 
compressed towards nuclear density. 

|Shternin et al.| ( |2007| l compared time-dependent calcula- 
tions of deep crustal heating and subsequent cooling to the 
observations of KS 1731-260. They found that a low value 
of thermal conductivity corresponding to an amorphous crust 
could be ruled out because it would give cooling on a much 
longer timescale than observed. In this paper, we confirm 
this result, and go further by describing the basic physics that 
sets the shape of the cooling lightcurve, and calculating in de- 
tail the constraints on Q\ mp and other crust parameters com- 
ing from the cooling lightcurves of both KS 1731-260 and 
MXB 1659-29. jCackett et al."1 ( [20Q6l > found that both of these 
decays could be fit with an exponential decay to a constant, 
although a single power-law (L oc r a , with a = 0.50 ± 0.03) 
also adequately fits the data for KS 1731-260 (Cackett et al. 



2008 ). We show here that the lightcurve of a cooling crust is 
a broken power law. The initial power law decay gives a di- 
rect measure of the temperature profile, and hence the thermal 
flux during outburst, in the outer crust. The time of the break, 
at hundreds of days post-outburst, corresponds to the thermal 
time where the solid transitions from a classical to quantum 
crystal, close to neutron drip. At late times, the luminosity 
levels off at a value set by the neutron star core temperature. 

We start in §2 by describing our time-dependent cooling 
calculations and an analytic model of the results, and go on 
in §3 to calculate the constraints on crust parameters coming 
from comparison with the observed cooling of KS 1731-260 
and MXB 1659-29. The Appendix discusses the details of 
our crust models. 

2. MODELS OF CRUST COOLING IN KS 1731-260 AND 
MXB 1659-29 

2.1. Hydrostatic structure of the crust 

Because the temperature is always low relative to the elec- 
tron and neutron Fermi energies, we can solve for the tem- 
perature and luminosity using a hydrostatic structure. In the 
crust, the pressure P makes a convenient Eulerian coordinate, 
and we integrate the equations ( |Thorne|1977| > for the radius r, 
gravitational mass M, and potential <p, 
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Here 1 + z = [1 - 2GM/(rc 2 )]- 1/2 , g = GM(l + z)/r 2 is the 
gravitational acceleration, and p is the density of mass-energy. 
We have neglected terms 0(pr i /Mc 2 ), as appropriate in the 
crust. As boundary conditions, we assume a transition den- 
sity to uniform npe matter at n — 0.08 firT 3 (consistent with 



recent studies of clustering in uniform nuclear matter; Oya- 
|matsu & Iida|2007] l, and set M and r a ccording to a neutron 
star model computed using the EOS of Akmal et al. ( 1998 ). 
We integrate outwards to a pressure P = 2.3 x 10 ZD ergs cm" , 
corresponding to a column depth from the surface 1 of Pfg = 
10 12 g cm 2 , at which point we apply the third boundary con- 
dition (f>(r = R) = (c 2 /2)ln[l - 2GM/(Rc 2 )]. The integration 

1 In a thin layer, the column depth is p dr as P/g; in this paper we will 
use the term to refer to y = P/g. 



is performed using a standard fourth-order Runge-Kutta al- 
gorithm, and the output is constrained to generate points uni- 
formly distributed in In P for use in the time-dependent code 
(§ 2.2 1. Our equation of state, as well as our model for the 
composition, is detailed in the Appendix. 
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Fig. 1. — Temperature in the neutron star ocean (T[,) a s a function of pho- 
tosphe re temperature r c (f {solid line). The relation of Gudmundss oiTet al.| 
1 1983 dashed l ine) is shown for comparison. We also show two models from 
Potekhin et al. 1 1997): their "fully accreted" model dotted line and a "par- 
tially accreted" model (dot-dashed line) in which the light elements are in the 
region P/g < W 9 gcm~ 2 . Note that for these latter two models, the temper- 
ature Tt is taken at a density 10 10 g cirT 3 (P/g a 4.3 X 10 13 g cm -2 ), which 
is somewhat deeper than the boundary used in our calculations. For compar- 
ison with these cases, we also show (thin solid line) our relation obtained by 
integrating to this depth. 



2.2. Time-dependent Heating and Cooling 

The time-dependent equations for the evolution of temper- 
ature and luminosity are 
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where e nuc and e v are the specific nuclear heating and neutrino 
emissivity, C is the specific heat, and K is the thermal conduc- 
tivity. We solve these equations using the method of lines. We 
use the common technique of defining Le 2 ^ at the midpoints 
of our grid by interpolating \nr 2 Ke^^ € ' /(l + z) and differenc- 
ing 

Te 4>/c- 

; as a result the divergence term in equation (pb is 
second-order and explicitly conserves flux. This procedure 
yields a set of coupled ordinary differential equations, which 
we then integrate using a semi-implicit extrapolation method 
(see |Press et al.| 1 992 and references therein). Our calculation 
of C, K, e nuc , and e v is described in the Appendix. 

We used two different boundary conditions for the core. 
The first is to simply assume a constant temperature, which 
we fit to observations. The second is to match the inwards 
luminosity at the crust-core boundary to the neutrino emis- 
sion from the core using a tabulated T c -L v relation for differ- 
ent assumptions of the core neutrino emissivity. In this way, 
we self-consistently solve for the core temperature appropri- 
ate for the assumed core physics rather than treat it as a free 
parameter. Unless the quiescent interval is long, we find that 
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TABLE 1 

Parameters for numerical integrations 



quantity 


unit 


KS 1731-260 


MXB 1659-29 


M 


gs~' 


10 17 


10 17 


T 00 
core 


10 7 K 


4.6 


2.6 


outburst duration 


yr 


12.0 


2.5 


recurrence time 


yr 


100 


21 



the core temperature is essentially constant over an outburst- 
quiescence cycle. 

The boundary condition at the surface is more ambiguous. 
During an outburst, the temperature in the neutron star en- 
velope is set by the burning of hydrogen and helium, and 
(possibly) fusion of light elements such as 12 C. Our code 
does not track this burning, and so we fix the temperature at 
P/g = 10 12 gcirT 2 at a fixed value Tf,. This column is roughly 
where superburst ignition occurs, and should demarcate the 
bottom of the region containing light element, unstable reac- 
tions. Because our first goal is to fit the lightcurve, for now we 
do not set T/, to any expected value a priori, but rather leave 
it as an adjustable parameter. 

During quiescence, we calculate the cooling flux at the top 
of the grid using a tabulated relation between and the 
temperature obtained by integrating the steady-state thermal 
structure of the neutron star envelope ( |Brown et al.|2002) . In 
these integrations, we fix the envelope to be pure 4 He down to 
a depth P/g = 10 9 g cirT 2 , with a layer of pure 56 Fe down to a 
depth P/g = 10 12 gem -2 . The resulting relation (Fig. |T| solid 
line) resembles that of Gudmunds son et aL] ( |1983| dashed 
li ne) at low T e g, but tren ds towards the fully accreted model 
of Potekhin et al. ( 1997 dotted line) at higher T^g. We also 
show the calculation of ( Potekhin et al.|199~7 1 for a light ele- 
ment envelope that is a accreted to a density (dot-dashed line). 
Note that our boundary is not quite as deep as the one used 
by |Gudmundsson et al."] ( |1983| ) or |Potekhin et al.| ( |1997) . For 
comparison, we therefore show (thin black line) our relation 
Th(T e ff) when the bottom of the 56 Fe layer is ta ken at th e same 
density as used by [Gudmundsson et al. ( 1983 1 and Potekhin 
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2.3. Numerical Results 

We now describe our method for fitting the data, using 
MXB 1659-29 as an example. With the hydrostatic struc- 
ture constructed (eqs. |2)-||4)), we use the surface gravity and 
redshift to tabulate Ti,(T^), where T/, is the temperature at 
y - Jtop- We set the core temperature from the last obser- 
vation of MXB 1659-29 (|Cackett et aTJ|2008l). Although 



our code can follow the thermal evolution of the core, the 
change in temperature is not large unless the quiescent inter- 
val is ~ 10 3 yr, the core cooling timescale. We therefore 
found it more convenient to fix the core temperature at this 
value. Table [T] shows the parameters for the calculations — 
accretion rate, core temperature, outburst duration, and recur- 
rence time — described in this section. In both cases, we used 
a neutron star with a mass 1 .62 M , radius 11.2 km, surface 
gravity 2.27 x 10 14 cms -2 , and surface redshift 1.32. 

We base the outburst and recurrence times on observa- 
tions. KS 1731-260 was first discovered in August 1989 with 
MzV/Kvant (Sunyaev 1989 ) and went into quiescence in 2001 
(Wijn ands et al.||20Ul) . Its recurrence time is unknown, al- 
though it has been quiescent since then. So long as the recur- 
rence time is longer than the thermal relaxation time of the 



crust (as determined by Shternin et al. 2007 and confirmed 
here), our results are not sensitive to the recurrence time since 
we are using the core temperature as a parameter in our fit. 
MXB 1659-29 has had two outbursts since being discovered 
in 1976 ( |Lewin et al.||1976] l. It was initially detected with 
SAS3 and HEAO through 1978; subsequent observations with 
a variety of instruments (see Wijn ands et al.|2003] and refer- 
ences therein) failed to detect the source until it was detected 
in outburst again in 1999 ( |in 't Zand et al. 1999[ l. This out- 
burst was monitored with RXTE/ASM until the source went 
into quiescence in September 2001. We adopt an outburst du- 
ration of 2.5 yr and a recurrence time of 21 yr. For both 
sources we adopt an accretion rate M = 10 17 gs" 1 , which 
is consistent with recent estimates (Gallowa y et al.|2008 [ see 



§ 3.2 1, for example from the X-ray luminosity using the re- 
lation for energy radiated by radially infalling material (see, 



Ayasli & Joss|l982| ) M = (1 + z)L x /c 1 z. There is un- 



e. g.,_ 

certainty in deriving the mass accretion rate from the source 
luminosity, and we shall explore the sensitivi ty of our conclu- 
sions to the assumed mass accretion rate in §|3.2| 
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Fig. 2. — An example of a lightcurve that fits the observed cooling of 
MXB 1659-29. The numerical model is shown as a solid curve and has 
gimp = 4.0 and T b = 3.8 X 10 s K. The dotted curve shows the correspond- 
ing toy model lightcurve from j |2.4| and we also show (grey dashed line) the 
slope given by eq. \\2\ . 

To generate a starting model, we first compute the thermal 
structure for a steadily accreting neutron star at the outburst 
accretion rate M. We then turn on the time-dependent terms in 
equations (BJ-([5J. We run through several outburst/quiescent 
cycles, and then generate a quiescent lightcurve with finer 
time resolution. Because the recurrence time is unknown for 
KS 173 1 -260, we use in the code an arbitrary recurrence time 
of 100 yr; note that our results do not depend on the value of 
the recurrence time, however, since the core temperature is 
held fixed. The code stores snapshots of the temperature dur- 
ing this "high resolution" run so that the physical properties 
of the crust during its gloaming can be reconstructed. 

Several runs with different values of Q lmp and 7j were then 
made, adjusting the chosen values to provide a good fit by eye 
to the lightcurve. Although there is no physical reason why 
Qi mp should be the same value throughout the crust, there are 
no reliable calculations of how gimp would evolve with depth, 
and we therefore choose to set Q lmp to a constant. Indeed, 
we shall show (§ 2.4 1 that it is really the value of gimp in the 



inner crust that is important for determining the lightcurve. 
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Fig. 3. — A model lightcurve for KS 1731-260, with Qj mp = 1.5 and 
T b = 2.5 X 10 8 K. 

Figure El shows a fit for MXB 1659-29, with gimp = 4.0 and 
Tj, — 3.8 x 10 8 K. The lightcurve is a broken power law with 
a break at t » 300 d, going to a constant at late times. 

Repeating this procedure for KS 1731-260, we find an ac- 
ceptable fit with 2imp = 1-5, Tb = 2.5 x 10 8 K, as shown 
in Fig. [3| No te that the errors on the temperature given by 
Cackett et al. ( 2006) are (erroneously) s tated to be 1-cr errors , 
but are in fact 90% confidence limits (Cac kett et al.] 2008 ). 
We show 1-cr errors in all figures in this paper. For the ob- 
servations of MXB 1659-29 we use the published 1-cr errors 
( |Cackett et aT1 [2008 ); for the observations of KS 1731-260 
we assume a Gaussian distribution to adjust the published 
90% confidence limits (Cackett et al. 2006) to 1-cr errors. 

2.4. Simple understanding of the lightcurve 

We now describe a simple model that accurately reproduces 
the lightcurve from the time-dependent calculation, and re- 
veals the basic physics underlying the lightcurve. This will 
allow us to understand the effect of different parameters on the 
lightcurve. A similar approach has been applied to lightcurves 
of white dwarfs cooling following a dwarf nova ( Piro, Ar- 
ras, & Bildsten 2005), and the early phase of superburst 
lightcurves (Cumming et al. 2006, Appendix A). 

We first note that during the long outbursts of 
MXB 1659-29 and KS 1731-260, the crust tempera- 
ture profile is very close to the thermal steady-state profile 
at the outburst accretion rate. This is shown in Fig. [4] for 
MXB 1659-29, in which we compare the temperature at 
the end of the outburst of MXB 1659-29 in our numerical 
calculations (top panel, dotted line) with the temperature 
profile of a steady-state calculation at the outburst accretion 
rate (top panel, dashed line). The largest difference (Fig. |4j 
bottom panel), in the inner crust where the strongest heat 
sources lie, is only 4% percent. Therefore a good approx- 
imation to the initial temperature profile for the cooling is 
the steady-state profile. This is an even better approxima- 
tion for KS 1731-260, which had a longer outburst than 
MXB 1659-29. Our calculation is therefore different from 



that of Ushomirsky & Rutledge (2001), who injected the 
entire heat deposition of the outburst at once. 

Starting with the initial temperature profile, we can un- 
derstand the evolution of the cooling layer and the resulting 
lightcurve by noting that at a given depth, the thermal evo- 
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Fig. 4. — (Top panel) Temperature in the crust of MXB 1659-29 at the end 
of a 2.5 yr outburst (dotted line) and that of a neutron star in thermal steady- 
state accreting at the outburst accretion rate (dashed line). Note that the core 
temperature and the temperature at a column y to p = 10 12 gem -2 are fixed. 
(Bottom panel) Ratio of the steady-state temperature to that at the outburst 
end. 

lution occurs on the characteristic thermal timescale associ- 
ated with that depth. This is illustrated in the middle panel 
of Figure[5] which shows snapshots of the temperature profile 
of the crust as it cools during quiescence. At a given time, 
the temperature profile has two parts: the inner layers have 
not yet started to cool and still have the temperature profile 
corresponding to the initial condition (the steady state profile 
during outburst); the outer layers have relaxed thermally and 
the temperature profile there corresponds to a constant out- 
wards flux. The transition occurs at a depth where the thermal 
time at that depth is equal to the current time. In the bottom 
panel of Figure [5] we show the thermal time as a function of 
depth, whe re we calculate the therm al time from the surface 
following Henyey & L'Ecuyer ( 1969 1, 



(7) 



where p is the density, Cp the specific heat, and K the thermal 
conductivity. In the top panel of Figure [5] we show the tem- 
perature profiles as a function of the thermal time. This shows 
directly that the deviation of each of the dashed temperature 
profiles away from the initial temperature profile occurs at a 
depth where the thermal time is approximately equal to the 
time since cooling began. 

The temperature of the inner crust is also affected by con- 
duction of heat into the core. We show the timescale for 
thermal diffusion into the core in Figure [5] (bottom panel, 
dotted line). The two thermal times intersect at a depth 
(P/g ~ 10 16 gem -2 ) where the thermal diffusion time is 
ss 400 d. After this point the temperature in the inner crust 
drops markedly (top panel). 

This understanding suggests a simple model of the 
lightcurve. We start with the initial temperature profile set by 
the steady-state profile at the outburst accretion rate. Then, for 
each time t, we locate the depth at which t = t . We then find 
the outwards flux in a constant flux solution that has a tem- 
perature equal to the initial temperature at the depth where 
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parameter T at which the ions crystallize (see Appendix A.l I. 
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Fig. 5. — Illustration of the cooling behavior of our analytical model. The 
bottom panel shows the thermal diffusion time, as a function of P/g, to the 
surface (eq. |9j; solid line) and to the core (grey dashed line), as well as the 
depth from the surface (dotted black line, right axis). The middle panel shows 
the temperature as a function of P/g. In the top panel, we plot the temper- 
ature (solid line) against t at the start of the quiescent phase. Subsequent 
"snapshots" of the temperature in the crust (grey dotted lines) at t = 3.1, 10, 
31, 107, 305, and 504 d plotted against r (top panel) and P/g (middle panel). 

t — t. This value of flux is the flux emerging from the surface 
at time t. The dotted curve in Figure [2] shows a lightcurve cal- 
culated in this way, using the same parameters (2imp, Tj,, and 
T c as the numerical model. The simple model shows excellent 
agreement with the numerical model. 

The origin of the broken power law nature of the lightcurve 
lies in the change in slope of the thermal time with depth that 
occurs close to neutron drip (see the lower panel of Fig. [5j 
neutron drip occurs at P/g w 5 x 10 15 g cm ). The decrease 
in slope is primarily due to the suppression of C p in the inner 
crust, shown in Figure [6] The ion contribution to the specific 
heat (bottom panel, dotted line) decreases on going to higher 
densities roughly as (T/&d) 3 , where the Debye temperature 

&d ® P - (h/ks) [47rZ 2 e 2 n; on /(Am u )J , the plasma tem- 
perature of the ions. We assume that the neutrons in the inner 
crust are superfiuid, in which case they have a negligible con- 
tribution to the heat capacity (Fig. [6] bottom panel, dot-dashed 
line). The thermal time also depends on the thermal conduc- 
tivity, which changes from being set by phonon scattering in 
the outer crust to impurity scattering in the inner crust {top 
panel, dashed line). Electron-electron scattering, although in- 
cluded in our calculations, is not a sign ificant com ponent of 
the total thermal conductivity (Sht ernin et al.|20"07) , and we 
do not show it in Fig.]6] The slight step in the ion specific heat 
at P/g < 10 13 gcm _2 (Fig.|6] bottom panel) is caused by the 
liquid-solid transition in the crust. Our code does not follow 
the crystallization front and hence does not include the latent 
heat. The depth where crystallization occurs is so shallow, 
however, that this omission does not appreciably affect the 
lightcurve, unlike the case for the cooling of white dwarf stars 
(see |Chabrier| 1 9 99 ; Hansen 2004, and references therein). We 
note that observations taken shortly (< 10 d) after the end 
of the outburst could potentially detect the effect of the latent 
heat; this would provide an independent constraint on the tem- 
perature in the crust and a check on the value of the plasma 
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Fig. 6. — Physical quantities setting the thermal diffusion time (cf. eq. (5)). 
The top panel shows the total thermal conductivity (dotted line), and the con- 
ductivity resulting from electron-phonon (solid line) and electron-impurity 
(dashed line) scattering. These quantities are computed for the temperature at 
the end of the outburst for the run corresponding to MXB 1659-29. The bot- 
tom panel displays the specific heat (solid line), in units of IcbNa/A, and the 
contributions from ions (dotted line), electrons (dashed line), and neutrons 
(dot-dashed line). The slight step in Cp (bottom panel) at P/g < 10 13 gcnT 2 
is where the ions crystallize. The top axis of the plot indicates the density as a 
function of P/g. The scale is not linear because of the change in the effective 
polytropic index in the inner crust. 

With some approximations, the same arguments allow us to 
make an analytic approximation to the lightcurve. The slope 
of the cooling curve can be written 



dlnT 



eff 



dlnf 



dlnr eff\/dlnr\/dlny 



dlnJ / \ dlny / ydlnr 



(8) 



The first factor on the right hand side is the slope of the r en -- 
T relation, d In T™ /din T * 0.45-0.63 (Fig.Q. The second 
factor is the temperature gradient in the initial model. The 
third factor is the dependence of thermal time with column 
depth. We can obtain this analytically by noting that dur- 
ing the early part of the lightcurve, when the cooling wave 
is in the outer crust, we can approximate the heat capacity as 
Cp ~ 3ks/(Am ll ) the classical heat capacity of a lattice, and 
use an a ppro x imat e expression for the phonon conductivity 
(see eq. | A3 |- |A4| } lj Inserting these expressions into the ex- 
pression for t (eq. |7|), we find 2 

5/4 



34 d-^ 4 



£)(£)(?)• » 



where gi4 = g/10 cms , giving d\nr/d\ny = 3/4. 

A measurement of the slope during the early part of the 
lightcurve directly measures, therefore, the temperature gra- 

2 In principle, we could include the dependence of Y e on y (eq. |A2| ), but 
for simplicity we shall leave it free in this formula. This expression agrees 
well with our numerical calculations (Fig. [5] bottom panel). 



(10) 



dient in the outer crust at the end of the outburst, 

dlnT ^ 15 din 7^. 

dlny 11 dlnf 

The fact that the effective temperature is decreasing with time, 
implies that temperature decreases with depth in the crust: the 
temperature profile in the outer crust is inverted during out- 
burst. If the temperature increased with depth in the crust, 
equation ( 10 1 implies that we would see an increasing temper- 
ature with time, and this can in fact be seen in the numerical 



models of |Rutledge et al. ( |2002| see their Fig. 3). 

A different way to write equation ( [T0| is in terms of the 
inwards flux in the outer crust, since the flux determines the 
temperature gradient. We write the flux as F — -KdT/dr and 
use the expression for K, eqs. (A3 i-( A4i to obtain 



dlnT 
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giving 
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where we have suppressed the other factors from equa- 
tion (fTTji and we use dlnT^/dlnT = 0.55. As an example, 
we plot a line representing this slope in Fig. [2] {grey dashed 
line), but using dlnT^/dlnT = 0.45, as appropriate at this 
r e ff. The slope of the early-time power law directly measures 
the flux in the rust. 
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Fig. 7. — Our fit to the lightcurve for KS 1731-260 (dotted line) compared 
to calculations for which the outer boundary condition during outburst was re- 
placed with the same 7i(T c ff) relation used during quiescence (dashed, solid 
lines). For the case where other parameters were held fixed (dashed line), 
the lightcurve falls markedly below the observed values. When the outburst 
accretion rate is increased to 3.5 X 10 17 g s -1 and the impurity parameter de- 
creased to 2imp = (solid line) a better agreement is found with the observed 
lightcurve. A clear distinction between these scenarios could be observed 
during the first fortnight of quiescence. 



Our inference that the temperature gradient in the crust is 
inverted, i. e., the temperature decreases with depth, also im- 
plies that an inward-directed heat flux enters the crust from 
the top. To show how this is critical for the early decay, Fig- 
ure U\ compares our original fit (dotted line) for KS 1731-260 
witn runs in which we do not hold the temperature at the 
top of the crust fixed, but rather use the same boundary 
condition (see Fig. [T| during both outburst and quiescence 



(dashed and solid lines). We first changed only the condi- 
tion on Ti„ while holding all other parameters fixed (dashed 
line). Although the break in the lightcurve occurs on the 
same timescale, the calculation with deep heating only is too 
faint to match observations. We then increased the accretion 
rate to M = 3.5 x 10 17 gs" 1 , which gives a crust heating 



rate comparable to that used by Shternin et al. (2007 1. Be- 
cause the higher temperature decreases the thermal conduc- 
tivity, we set Qi mp — to compensate. The resulting solution 
(solid line) gives a better fit, although its early time behavior 
still undershoots the first observation. For both latter cases 
(dashed and solid lines), the temperature profile is inverted 
for P/g > 10 13 gem" 2 ; that is, most of the heat produced in 
the outer crust during outburst is conducted inward. As a re- 
sult, is decreasing with time starting a few days after the 
end of outburst in these cases. 

As noted by Shternin et al.| (|2"007), the observed lightcurve 
of KS 1731-260 can be explained without recourse to an 
inward-directed heat flux provided that M is larger than the 
value we assume in this paper (cf. Fig. [7}. We have made the 
same comparison for MXB 1659-29 (Fig. |8j: we used the 
same Tb(T e g) relation during both outburst and quiescence, 
and then increased M while decreasing Qimp- In addition to 
our best-fit solution with fixed Tb (dotted line), we show the 



case for M 
M = 9 x 10 17 gs 



5 x 10 17 gs" 1 



and Q 



imp 



1 (dashed line) and 



and gimp = 0. Without a fixed Tb during 



outburst (or equivalently a sizeable inward-directed heat flux 
in the shallow outer crust), we do not find an acceptable nu- 
merical solution consistent with the observed MXB 1659-29 
lightcurve. 
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Fig. 8. — A comparison between the observed lightcurve of MXB 1659-29 
and several trial numerical solutions: 1) our best-fit solution (dotted line; 
cf. Fig. [2) with T b = 3.8 X 10 8 K and M = 10 17 gs -1 ; a solution with 
M = 5 X 10 17 gs~\ Qtap = 1, and T b = T b (T zS ) (dashed line); and a 
solution with M = 9 X 10 17 gs -1 , Qimp = 0, and T b = T b (T sS ). 



To illustrate this point, we computed, using equation ( 12 1, 
the required flux necessary to match the observed power-law 
slope between the first two observations. For KS 1731-260 
the first two observations imply a flux of 0.7 MeV/;« u for 
M = 10 17 gs 1 . For MXB 1659-29, the flux required 
to match the observed power-law slope is 0.8 MeV/m u for 
M = 10 17 gs" 1 . Our numerical solutions are consistent 
with this estimate: our numerical model of the KS 1731-260 
lightcurve (Fig. [3j has a local inward flux of 0.5 MeV/m u , 
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while our model of the MXB 1659-29 lightcurve (Fig. [2} has 
1.1 MeV/m u . This value of the flux is well above the es - 
timates available from electron captures ( Gupta et al.||2007] l. 
Moreover, this flux would have to emanate from a depth for 
which the thermal time is less than the time of the first ob- 
servation, which is Pig < 2 x 10 14 gem -2 , corresponding to 
p < 3 x 10 10 g cirT 3 . This is a lower density than that of many 



crust electron captures (Gupta et al. 2007 Haensel & Zdunik 
2008) > and e ven that of light elemen t pycnonuclear reactions, 
such as 24 ( |Horowitz et al.|2008b] ). Hence, although a larger 
M reduces the total amount of heating per nucleon required, 
matching the first data point is still difficult. Observations of 
both sources are too sparse on timescales < 10 2 d, however, to 
draw firmer conclusions about the existence and nature of this 
heating. Observations occurring within the first two weeks af- 
ter the outburst ends are critical for constraining the depth and 
strength of heat sources in the outer crust. 

3. CONSTRAINTS ON THE MODEL PARAMETERS 

In Figures |2]and[3] we show models that fit the data well for 
MXB 1659-29 and KS 1731-260. We now address the un- 
certainties in the fitted model parameters. The cooling mod- 
els have a range of physics input. To allow us to investigate 
the full range of parameter space, we adopt the simplified 
model of the cooling curves described in §2.4 as it allows us 
to rapidly generate a lightcurve without running the full time- 
dependent simulation. 

To calculate the constraints on model parameters, we adopt 
a Markov Chain Monte Carlo (MCMC) method, using the 
Metropolis-Hastings algorithm to generate a sequence of sam- 
ples from the posterior probability distributions (e. g., |Gre- 
|gory|20 05 ). Although we have a small number of parameters 
which would allow a grid search of the chi-squared space, 
we prefer the MCMC method for its simplicity in implemen- 
tation and generating marginalized probability distributions. 
We run the chains multiple times from different starting points 
to check the convergence and robustness of the resulting prob- 
ability distributions. Typically, we find that ~ 10 4 samples 
are adequate, although we have run the chains longer to check 
convergence. 

Although the simplified model lightcurves agree well with 
the numerical results (as shown for a specific example in Fig- 
ure [2}, the agreement is not exact, and as a result the best-fit 
parameters derived from the two models are different. For 
example, whereas the best fit lightcurve, for MXB 1659-29, 
from the time-dependent simulation has Qi mp = 4 (Fig.|2j), the 
simplified model gives a best fit value Q\ mp = 2. In this sec- 
tion our focus is, however, on how well we can constrain the 
main parameters of the model, and how each of these param- 
eters changes the lightcurve, rather than the best-fit values. 
By running a sample of comparison models, we have checked 
that the uncertainties in the parameters derived from the sim- 
plified model and time-dependent simulations are similar. 

3.1. Constraints for a given hydrostatic structure 

For a given hydrostatic structure, which is set by the neu- 
tron star mass and radius, the three main parameters that affect 
the lightcurves are the temperatures at the top of the crust 7j, 
the core temperature T c , and the impurity parameter Qimp- We 
assume for now that the accretion rate, which determines the 
overall heating rate in the crust (see Appendix), is determined 
from the observed X-ray luminosity. We fix the neutron star 
mass M = 1.6 M Q and radius R = 11.2 km, and the outburst 
accretion rate M = 10 17 g s _1 for both sources. These values 



match the parameters used in the time-dependent code to cal- 
culate the models shown in Figures [2] and [5] The correspond- 
ing redshift factor is 1 + z — 1.32 and gravity is g\n = 2.3. 

Figure [9] shows the resulting probability distributions for 
Gimp, and Tt,. The values of T™ and Tb are well- 
constrained in each case. These two temperatures are con- 
strained by the need to match the first and last data points in 
each set of observations. The first data point gives a snapshot 
of the temperature near the top of the crust (thermal time of 
tens of days), whereas the last point gives a measurement of 
the core temperature (we find that the lightcurve has leveled 
off in both sources in our models). 




Fig. 9.— Probability densities for g imp , Tf and T b for MXB 1659-29 and 
KS 1731-260. The peak of the probability distribution is marked by a cross, 
and the contours show the 68% and 95% confidence intervals. The prior 



probability for gimp is taken to be uniform in log from Q n 



10" 3 to 10 2 . 



For each source, we fix the outburst accretion rate at M = 10 g s and the 
neutron star mass and radius at M = 1.6 M Q and R = 11.2 km. 




10 d 10 2 
time since outburst end (d) 

Fig. 10. — Lightcurves of MXB 1659-29, for different choices of the im- 
purity parameter Qimp in the crust. In both case, we show the best fit model 
(solid line). The other solutions have Q Imp = (dot-dashed line), 1 (dotted 
line), and 10 (dashed line). The left-hand panel shows the case for which all 
other parameters are held constant; in the right-hand panel the temperature at 
a column y lo „ = 10 12 gcnT 3 was adjusted so that all solutions matched the 
first data point. 



The impurity parameter £}; mp is very tightly constrained in 
MXB 1659-29. To understand why this is the case, we show 



in Figure 10 several lightcurves from our time-dependent sim- 



ulations fofMXB 1659-29 with different values of 2imp- In 
the left panel, we keep all other parameters fixed as we vary 
2imp; in the right panel, we adjust Tb to match the first data 
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point as we vary Qi mp . The effect of increasing Q lmp is to de- 
lay the cooling. This can be understood in terms of the initial 
temperature profile at the end of the outburst. When the crust 
has a larger impurity level, the inner crust must be hotter to 
be able to conduct the heat released in the crust into the core. 
This increase in the inner crust temperature leads to a hotter 
outer crust, with a shallower temperature gradient, giving a 
lightcurve that falls less quickly. 

As Figure [9] shows, the correlations between the fitted pa- 
rameters are such that larger Qimp values lead to lower values 
of T c and 7j, i.e. to compensate for the delayed cooling due to 
increase in Qimp, the overall temperature scale set by T c and 
Tj, decreases. In KS 1731-260, the probability distribution of 
(2imp has a peak at a similar value to MXB 1659-29, but with 
a long tail to small values of Qi mp - In fact, as can be seen in 
Figure 1 1 the fits are not sensitive to the impurity parameter 



for (2i mp < 1, which results in a flat probability distribution in 
log Q\ mp , reflecting the assumed prior. For both sources, Q; mp 
values larger than 10 are ruled out. 

For MXB 1659-29, we have used the temperatures derived 
by Cackett et al. (2008) assuming a distance to the source of 
10 kpc. In that paper, spectral models for different distances 
d — 5 and 13 kpc are considered, which leads to a system- 
atic decrease or increase in the effective temperatures by 10- 
20%. The reason that the fitted effective temperatures depend 
on distance is that the peak of the thermal spectrum lies out- 
side the X-ray band, making the fitted temperature sensitive to 
the overall luminosity scale. To investigate the effect of such 
systematic variations, we have calculated the constraints on 
the models with the effective temperatures for MXB 1659-29 
all decreased or increased by 20%. The effect is to change 
the central value of each distribution by up to 50%, with the 
width staying about the same. The conclusion that Qimp is of 
order unity is unaffected by these systematic variations. 




Fig. 11. — The joint probability distribution of Q\ mp and M for 
MXB 1659-29 (dotted contours) and KS 173 1-260 (solid contours). In each 
case, the peak of the probability distribution is indicated by a cross; the two 
contours enclose 68% and 95% of the probability. 



Appendix for details). The calculations so far have taken a 
fixed accretion rate M = 10 17 g s _1 . Instead, we now calcu- 
late the constraints on M assuming a uniform prior probabil- 
ity for M between (i.e. no deep heating) and 10 18 g s _1 (ten 
times our fiducial rate). The results are shown in Figure[TT] in 
which we give the derived joint probability distribution for M 
and Q imp for each source. The temperatures Tf, and T c are not 
sensitive to variations in M, since they are essentially fixed by 
the first and last observed values of T°^, 

For both sources, we find an anti-correlation between M 
and Q; mp in the best-fitting solutions. The explanation for the 
anti-correlation is that an increased M gives an increased heat- 
ing rate, making the inner crust hotter. To compensate for this, 
Qimp must decrease, cooling the inner crust by making it eas- 
ier for heat to be conducted into the core. 

The values of M derived from the cooling curves match 
well with the accretion rates derived from observations 
of the persistent X-ray luminosity during outburst. For 
MXB 1659-29, the range of flux observed during the out 
burst was * (0.4-1) x 10~ 9 ergs s" 1 cnr 2 (2.5-25 keV) ((Gal- 



|loway et al.|2008) . |Galloway efaiT| ( p008] l found a distanceof 
12 + 3 kpc for this source, assuming that the peak flux of pho- 
tospheric radius expansion bursts corresponds to the pure he- 
lium Eddington luminosity. Taking this distance and assum- 
ing a bolometric correction of a factor of 2, typical for these 
sources, gives M w (0.7-1.8) x 10 17 erg s _1 . The agreement 
with the constraints from the cooling curve is good, although 
lower than the maximum of the probability distribution for M. 

For KS 1731-260, [Galloway et al7| ( [2008| give a range of 
bolometric flux 1.6-10 x 10~ v erg cm s~\ which for their 
distance 7.2 + 1 kpc gives a range of accretion rates during 

A separate check on this 
a 

very regular sequence of X-ray bursts was seen, similar to 
the source GS 1826-24, which is known for being a very reg- 
ular burster. Assuming an ignition mass of 10 21 g for these 
regular bursts, which had a recurrence time of 2.59 + 0.06 h, 
we find M = 1.1 x 10 17 g s _1 , consistent with the X-ray flux. 
During the final > 1 year of the outburst, the flux was in the 
lower end of the flux range quoted earlier, so we expect the 
relevant value for the crust heating at the end of the outburst 



outburst of 0.5-3 x 10 17 g s" 1 
value is that at a flux level of 2.1 x 10~ 9 erg cirT 2 s _1 
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to be < 10 17 g s _1 , in good agreement with Figure 

An interesting aspect of our results is that both sources al- 
low solutions with low accretion rates much smaller than the 
accretion rates derived from the X-ray observations. Assum- 
ing that the observed accretion rate is within a factor of two 
of the true accretion rate onto the neutron star, this means that 
both cooling curves are consistent with a much lower amount 
of deep crustal heating than assumed in our models. In these 
models, however, a lower level of deep crustal heating from 
reactions in the crust is compensated by a larger inward heat 
flux from the neutron star ocean, because Tt, is held fixed. In 
reality, the physics of the implied, unspecified heat source in 
the neutron star ocean that supplies this flux also depends on 
the accretion rate. 



3.2. The accretion rate or overall heating rate in the crust 

The accretion rate M sets the overall amount of heating in 
the crust during the outburst. There are uncertainties in de- 
riving M from the observed X-ray luminosity, and in addi- 
tion, the amount of heating in the crust may differ from the 
1 .7 MeV per nucleon that we assume in our calculation (see 



3.3. Effect of varying neutron star mass and radius 

We now include the uncertainty in all six parameters of the 
model. The prior probabilities for T™ , Tt, Qi mp , and M are as 
described previously. In addition, we assume a uniform prior 
in the range 8-16 km for R and 1.1 to 2.5 M Q for M . The fit 

1.4 M and R = 10 km ( |Cackett 



for T™ used a fixed M 



et al. 2006 2008) because the quality of the data is insum- 
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Fig. 12. — Constraints on r™, 7i, and Qi rap , as in Figure|9] 
ing variations in M, R, and M. 



but now includ- 



cient to constrain these parameters independently; therefore, 
the constraints we derive should be considered indicative. The 
resulting constraints on T™, Tb, and Q lmp are shown in Fig- 
ure 12 This Figure should be compared with Figure [9] for 
which M, R, and M were assumed to be known. Including the 
uncertainty in the extra parameters broadens the probability 
distributions of T™, T/, and Q\ mp . The largest effect is on the 
probability distribution of Qi mp . For example, the constraint 
on (2imp is considerably relaxed for MXB 1659-29. For both 
sources, however, Q lmp values greater than several are ruled 
out even with the additional parameters included. The cen- 
tral values of T™ and Tb are similar to the values previously 
derived. 

The sensitivity of the derived value of (jimp on M and R 
is illustrated for MXB 1659-29 in Figure [T3] (we see the 
same effect for the KS 1731-260 data). We show the derived 
probability distribution for 2imp for three different choices 
of neutron star mass and radius. In each case, we keep the 
accretion rate fixed at our fiducial value M = 10 17 gs _1 . 
The allowed values of Q lmp increase with increasing surface 
gravity. This can be understood by considering the thermal 
time from a given density to the surface, which depends on 
the thickness of the layer and therefore varies with surface 
gravity (Lattimer et al. [1994). Rewriting the integral for the 
thermal time, equation (|9J, as an integral over pressure gives 
r°° oc (1 + z)/g 2 °c R 4 M~ 2 (1 + z) _1 . An increase in surface 
gravity shortens the cooling time, and Q lmp must increase to 
bring it back into agreement with the observed value. 

Figure 14 shows the joint probability density for M and 
R for each source. Although M and R are only weakly 
constrained, we see that the best-fitting values of M and R 
are correlated. The mass and radius enter the calculation 
of the lightcurve in several places besides the thermal time 
r°°. The relation between crust temperature and T°% de- 



pends on the surface gravity; for a fixed crust temperature, 
-^eff K 8 l ^/(l The initial temperature profile also changes 
with gravity. Using the Newtonian equations for the steady- 
state thermal profile, we see thatdr/dP = (l/g)(3KF/4acT 3 ), 
dF/dP = -e/g, so that the increase in flux due to the deep 
heating is smaller by a factor g, and the change in tempera- 
ture for a given flux is smaller by a factor g. The combina- 
tion of these different effects results in the observed correla- 
tion between the best fitting values of M and R. By inspec- 
tion we find that the slope of the relation is well-described by 
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Fig. 13.— The probability distribution of Q irap derived for MXB 1659-29, 
for three different choices of neutron star mass and radii. Left to right, in 
order of increasing surface gravity, they are (i) M = 1.4 Mo, R = 13 km, 
£14 = 1.4, l+z = 1.21 (ii)M = 1.6 M ,/? = 11.2 km, g M = 2.3, 1+z = 1.32 
and (iii) M = 2 M G , R = 10 km, g M = 4.2, 1 + z = 1.57. In each case, 
the accretion rate is fixed at our fiducial value M = 10 17 g s~'. Note that 
the spectral fits used to obtain T°l assume a fixed value of M = 1.4 M Q and 
R = 10 km. 




Fig. 14. — Constraint on the neutron star mass and radius. We assume a 
constant prior in mass between 1.1 and 2.5 M Q and in radius between 8 and 
16 km. The peak of the probability distribution is marked with a cross, and 
the contours enclose 68% and 95% of the probability. Note that the spectral 
fits used to obtain assume a fixed value of M = 1.4 M Q and R = 10 km. 



4. DISCUSSION AND CONCLUSIONS 

We have presented numerical simulations of the cool- 
ing of the neutron star crust in both KS 1731-260 and 
MXB 1659-29 following the end of long accretion outbursts. 
Our main results are as follows. 

1 . The lightcurve of a cooling crust is a broken power-law 
going to a constant at late times. The luminosity at late 
times is set by the neutron star core temperature. The 
slope of the early part of the lightcurve provides a direct 
measure of the flux in the outer crust during outburst 
(eq. (12)). The time of the break is set by the transition 
from a classical to quantum crystal, close to neutron 
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drip. The good fit of our models to the data suggests 
that the neutrons in the inner crust do not contribute sig- 
nificantly to the heat capacity, as expected if they were 
superfiuid. Observations of cooling quiescent neutron 
stars thus provide new evidence for the existence of a 
neutron superfiuid throughout the inner crust. 

As our models show, the observations to date probe 
the thermal relaxation timescale of the inner crust. 
The cooling timescale increases with increasing Qi mp , 
potentially giving a tight constraint on this parame- 
ter. The fits to the lightcurves of MXB 1659-29 
and KS 1731-260 both r equire Qj mp < 1 , in a gree- 
ment with the result of |Shternin et al.| ( |2007 1 for 
KS 1731-260. For our fiducial model, which has neu- 
tron star parameters M — 1 .6 M Q , R = 1 1 .2 km, and 
outburst accretion rate M = 10 17 g s _1 , the best fit val- 
ues are g imp = 4 for MXB 1659-29, and gimp = 1.5 for 
KS 1731-260. Reducing the surface gravity or increas- 
ing the accretion rate allows smaller values of Qi mp . 
Impurity scattering sets the thermal conductivity of the 
inner crust, and so our measurement of gimp refers to 
the conductivity in the inner crust, particularly close to 
neutron drip where the thermal time corresponds to the 
time of the break in the cooling curve. Interestingly, 
the values of Q lmv derived for both KS 1731-260 and 
MXB 1659-29 are very similar, and may indicate that 
a robust outcome of nuclear processing in an accreted 
crust is an inner crust impurity parameter of order unity, 
as suggested by calculations of nuclear transitions in the 
inner crust ( |Jones|2005||Gupta et al.|2008| ). 

The flux required to match the power-law slope be- 
tween the first and second observations is much larger, 
however, than expected from models of deep heating 
(Gupt a et al.l|2007| |Haensel & Zdunik|[20U8| . For 
KS 1731-260, we find that the lightcurve at > 100 d 
post-outburst can be fit using standard models of deep 
heating, if the accretion rate is larger than our esti- 
mate of 10 17 gs" 1 , in agreement with the findings of 
|Shternin et al. ( |2007| l. We do not find such a solution 



for MXB 1659-29, however; indeed setting the outer 
boundary condition to Tt, = Tb(T e g) drives the outburst 
accretion rate to roughly the Eddington limit. To obtain 
an adequate fit to the data, we require the temperature 
in the outer crust to be decreasing inwards, implying 
that an inward-directed heat flux enters the crust from 
the top. Moreover, our Markov Chain Monte Carlo fits 
with our simplified model (§ [3]l find best fit solutions 
with rather large values of T/,, so that the temperature 
profile is inverted. Our interpretation therefore differs 
slightly from that of |Shternin et al. (2007 ). This in- 
terpretation depends, however, on the first data points 
in each cooling curve, and so could be relaxed if these 
data points are contaminated, by residual accretion for 
example. Our calculations show that observations taken 
within the first two weeks following extended outbursts 
are ideal for mapping out the nuclear heating in the 
outer crust. It is this shallow heating that is most crit- 
ical for determin ing the ignition depth of superbursts 
(Gupta et al. 2007} . We shall investigate the heating re- 
quired to maintain the inferred high Tb, along with its 
implications for nuclear processes in the neutron star 
ocean, in a follow-up paper. 
We thank Chuck Horowitz, Ed Cackett, Nathalie Degenaar, 
Sanjay Reddy, Andrew Steiner, Lars Bildsten, Gil Holder, 
Bob Rutledge, and Sanjib Gupta for helpful discussions. 
EFB and AC acknowledge the hospitality of the Institute 
for Nuclear Theory, where this work took shape, and by 
the support of the Joint Institute for Nuclear Astrophysics 
(JINA) under NSF-PFC grant PHY 02-16783. EFB is sup- 
ported by the National Aeronautics and Space Administra- 
tion through Chandra Award Number TM7-8003X issued by 
the Chandra X-ray Observatory Center, which is operated by 
the Smithsonian Astrophysical Observatory for and on be- 
half of the National Aeronautics Space Administration un- 
der contract NAS8-03060, by NASA award NNX06AH79G, 
and by NASA ATFP grant NNX08 AG76G. AC acknowledges 
support from an NSERC Discovery Grant, FQRNT, and the 
Canadian Institute for Advanced Research (ClfAR), and as 
an Alfred P. Sloan Research Fellow. 



REFERENCES 



Aguilera, D. N., Cirigliano, V., Pons. J. A., Reddy, S., & Sharma, R. 2008, 
ArXiv e-prints 

Ainsworth, T. L., Wambach, J., & Pines, D. 1989, Phys. Lett. B, 222, 173 
Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Phys. Rev. C, 58, 
1804 

Altamirano, D., Wijnands, R., Degennar, N, Casella, P., Homan, J., Belloni, 

T., & Campana, S. 2007, The Astronomer's Telegram, 1178, 1 
Ayasli, S., & Joss, P. C. 1982, ApJ, 256, 637 

Bassa, C, Jonker, P. G., Nelemans, G., Steeghs, D., Torres, M. A. P., Kuiper, 
L., in't Zand, J. J. M., Rea, N., Maccarone, T., Kuulkers, E., Grindlay, J., 
Wijnands, R., & Mendez, M. 2008, The Astronomer's Telegram, 1575, 1 

Brown, E. F. 2000, ApJ, 531, 988 

Brown, E. F, Bildsten, L., & Chang, P. 2002, ApJ, 574, 920 

Cackett, E. M., Wijnands, R., Linares, M., Miller, J. M., Homan, J., & Lewin, 

W. H. G. 2006, MNRAS, 372, 479 
Cackett, E. M., Wijnands, R., Miller, J. M., Brown, E. R, & Degenaar, N. 

2008, ApJ, 687, L87 
Caillol, J. M. 1999, J. Chem. Phys., Ill, 6538 
Chabrier, G. 1993, ApJ, 414, 695 

Chabrier, G. 1999, in 11th European Workshop on White Dwarfs, ed. J.-E. 

Solheim & E. G. Meistas (San Francisco: Ast. Soc. Pacific), 369 
Chabrier, G., & Potekhin, A. Y. 1998, Phys. Rev. E, 58, 4941 
Degenaar, N., Wijnands, R., Wolff, M. T, Ray, P. S., Wood, K. S., Homan, 

J., Lewin, W. H. G., Jonker, P. G, Cackett, E. M., Miller, J. M., & Brown, 

E. F. 2009, MNRAS, in press 
Dewitt, H, & Slattery, W. 1999, Contributions to Plasma Physics, 39, 97 



Farouki, R., & Hamaguchi, S. 1993, Phys. Rev. E, 47, 4330 

Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 

2008, ApJS, 179, 360 
Gnedin, O. Y., Yakovlev, D. G., & Potekhin, A. Y. 2001, MNRAS, 324, 725 
Gregory, P. C. 2005, ApJ, 631, 1198 

Gudmundsson, E. H, Pethick, C. J., & Epstein, R. I. 1983, ApJ, 272, 286 
Gupta, S., Brown, E. F, Schatz, H., Moller, P., & Kratz, K. 2007, ApJ, 662, 
1188 

Gupta, S. S., Kawano, T, & Moller, P. 2008, Physical Review Letters, 101, 
231101 

Haensel, P., & Zdunik, J. L. 1990, A&A, 227, 431 

— . 2003, A&A, 404, L33 

— . 2008, A&A, 480, 459 

Hansen, B. 2004, Phys. Rep., 399, 1 

Henyey, L., & L'Ecuyer, J. L. 1969, ApJ, 156, 549 

Homan, J., van der Klis, M., Wijnands, R., Belloni, T, Fender, R., Klein- 
Wolt, M., Casella, P., Mendez, M., Gallo, E., Lewin, W. H. G., & Gehrels, 
N. 2007, ApJ, 656, 420 
Horowitz, C. J., Berry, D. K., & Brown, E. F. 2007, Phys. Rev. E, 75, 066101 
Horowitz, C. J., Caballero, O. L., & Berry, D. K. 2008a, Phys. Rev. E, in 
press 

Horowitz, C. J., Dussan, H, & Berry, D. K. 2008b, Phys. Rev. C, 77, 045807 
in 't Zand, J., Heise, J., Smith, M. J. S., Cocchi, M., Natalucci, L., & 

Celidonio, G. 1999, IAU Circ, 7138, 1 
Itoh, N, & Kohyama, Y. 1993, ApJ, 404, 268 
Jones, P. B. 2005, Phys. Rev. D, 72, 083006 



11 



Jonker, P. G., Bassa, C. G., Nelemans, G, Juett, A. M., Brown, E. E, & 

Chakrabarty, D. 2006, MNRAS, 368, 1803 
Lattimer, J. M., van Riper, K. A., Prakash, M., & Prakash, M. 1994, ApJ, 

425, 802 

Leinson, L. B., & Perez, A. 2006, Phys. Lett. B, 638, 1 14 
Levenfish, K. P., & Yakovlev, D. G. 1994, Astonomy Reports, 38, 247 
Lewin, W. H. G, Hoffman, J. A., Doty, J., & Liller, W. 1976, IAU Circ, 
2994, 2 

Oyamatsu, K., & Iida, K. 2007, Phys. Rev. C, 75, 015801 
Piro, A. L., Arras, P., & Bildsten, L. 2005, ApJ, 628, 401 
Potekhin, A. Y., Baiko, D. A., Haensel, P., & Yakovlev, D. G. 1999, A&A, 
346, 345 

Potekhin, A. Y, & Chabrier, G. 2000, Phys. Rev. E, 62, 8554 
Potekhin, A. Y, Chabrier, G, & Yakovlev, D. G. 1997, A&A, 323, 415 
Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, 

Numerical Recipes in FORTRAN (Cambridge: Cambridge Univerisity 

Press) 

Rutledge, R. E., Bildsten, L., Brown, E. E, Pavlov, G. G, Zavlin, V. E., & 

Ushomirsky, G. 2002, ApJ, 580, 413 
Sakano, M., Koyama, K., Murakami, H., Maeda, Y, & Yamauchi, S. 2002, 

ApJS, 138, 19 

Schatz, H., Bildsten, L., Cumming, A., & Wiescher, M. 1999, ApJ, 524, 1014 
Schinder, P. J., Schramm, D. N, Wiita, P. J., Margolis, S. H.. & Tubbs, D. L. 
1987, ApJ, 313, 531 



Sedrakian, A., Miither, H., & Schuck, P. 2007, Phys. Rev. C, 76, 055805 
Shternin, P. S., Yakovlev, D. G, Haensel, P., & Potekhin, A. Y. 2007, 

MNRAS, 382, L43 
Steiner, A. W., & Reddy, S. 2009, Phys. Rev. C, 79, 015802 
Sunyaev, R. 1989, IAU Circ, 4839, 1 
Thorne, K. S. 1977, ApJ, 212, 825 
Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 
Urpin, V. A., & Yakovlev, D. G. 1980, Soviet Ast, 24, 126 
Ushomirsky, G, & Rutledge, R. E. 2001, MNRAS, 325, 1 157 
Wijnands, R., Guainazzi, M., van der Klis, M., & Mendez, M. 2002, ApJ, 

573, L45 

Wijnands, R., Homan, J., Miller, J. M.. & Lewin, W. H. G. 2004, ApJ, 606, 
L61 

Wijnands, R., Miller, J. M., Groot, P. J., Markwardt, C, Lewin, W. H. G, & 

van der Klis, M. 2001, ApJ, 560, L159 
Wijnands, R., Nowak, M., Miller, J. M., Homan, J., Wachter, S., & Lewin, W. 

H. G. 2003, ApJ, 594, 952 
Yakovlev, D. G, Kaminker, A. D., Gnedin, O. Y, & Haensel, P. 2001, 

Phys. Rep., 354, 1 



APPENDIX 
CRUST MICROPHYSICS 

In this Appendix, we describe how we calculate the different pieces of microphysics that go into the crust models. 



Equation of state and composition profile 

We compute the equation of state for the degenerate, relativistic electrons by interpolation from tables of the free energy 
(Timmes & Swesty 20001 . For the ions we compute the EOS using the free energy fit of Cha brier & Potekhin] ( [19 98 ) for the 
liquid state and the fit of Farouki & Hamaguchi ( 1993| l for the crystalline state. With this choice of fits, the freezing point 
occurs where the free energies in the liquid and solid phase are equal, namely F = Z 2 e 2 /(akT) = 178, where a is the ion-sphere 
radius. Note that our fit to the liquid-phase free energy does not i nclude corrections (Potekhin & Chabrier 2000) based on recent 
Monte Carlo simulations (Dewitt & Slattery 199 9l|Caillol|1999|). As a result, our liquid-phase free-energy is slightly larger (by 
0.005% at F = 175) than the ex pression of |Potekhin & Cha brier (2000), and our freezing point is slightly higher than the value 
of F = 175 ± 0.4 determined by Potek hin" & Chabrier| ( |2000| l. Although we use the value of T at the freezing point determined 
by the ion free energy for a one-component plasma (OCP), it is im portant to note that the fre ezing point may differ significantly 
from this value because of polarization of the electron background (Potekhin & Chabrier 2000) and because the crust consists of 
a mixture of isotopes. Indeed, recent molecular dynamics simulations ( |Horowitz et al.|20Q7| > find that the freezing temperature 
for a mixture appropriate for X-ray burst ashes is smaller than the OCP, with T « 250 at freezing. 

We includ e the Debye re duction in the ion specific heat using an interpolation formula, and compute the Debye temperature 
according to Chabrierjl 1993| ). In the inner crust, we compu te the neutron specific heat as that of a degenerate gas, with suppression 
due to the pairing interaction ( [Levenfish & Yakovlev 1994| . For the 1 Sq pairing of free neutrons in the inner crust, we approximate 
the cri tical temperature as a Gaussian in the neutron Fermi wavevector, which approximates the calculation of Ains worth et al.| 
( 1989). The crust temperatures achieved in the models presented here lie below the critical temperature throughout most of the 
crust, and the neutrons do not contribute significantly to the total specific heat. 

The composition of accreting neutron star crusts has been calculated previously by Haensel & Zdunik (2008 ) and Gupta et al. 
(2007). They find that the composition of the crust changes abruptly with depth, at locations corresponding to thresholds for 



electron capture or pycnonuclear reactions. Rather than model the neutron star crust as a series of distinct layers, we instead fit 
the composition so that Z and A are smoothly varying functions of p. This approximation has a negligible effect on our results, 
but makes the integration simpler as it avoids jumps in the thermal properties. In the outer crust, we set the mass number A and 
compute Z in order to maintain /^-equilibrium. Using a simple liquid-drop model for the nuclear binding energy, 



B(A,Z) = a v A - a s A lp - a A 



(N - Z) 2 



A l/3' 



(Al) 



and minimizing the nucleon-specific Gibbs free energy G - Y e fj. e - B(A,Z)/A of a cell containing a nucleus and Z electrons with 
respect to Z for a fixed A, we obtain the electron fraction Y e in the outer crust as a function of electron chemical potential p e , 



\ 4a A j 



(A2) 



where we take 3 a A = 23.43 MeV and a c = 0.715 MeV. This formula accurately reproduces Y e in the crust models computed by 



Gupta et al. (2007). 



3 We obtain values for and ac by fitting experimental binding energies using the calculator at 



http : // 12 8 . 95 . 9 5 . 6 l/~intuser/ld . html 
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In the inner crust, the relation between chemical potentia l, Y e , and the abundance of free neutrons Fjy becomes more compli- 
cated. We find that the tables in Haensel & Zdunik ( 2008[ ) are adequately fit by defining A t o t as the total number of nucleons 



(including free neutrons) per nucleus and setting Z/A tot oc p~ 2/3 and A tot increasing as A tot <x /? 2/3 up to a maximum value set by the 
total number of pycnonuclear reactions that occur in the crust. We set a floor to Z/A tot of 0.03, and take A = max(A outel -, 0. 15A tot ), 
where A outer is the mean mass number in the outer crust. Figure [15] shows an example of a crust composition with A outer = 5 6 and 

A - 

^Mot.max — 

circles). 



896 in the inner crust, where we choose A tot max to agree with the calculation of |Haensel & Zd unik ( 2008 , Fig. [T3] open 
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Fig. 15. — Model of composition in the crust: nuclear charge (top p anel, dashed line) and mas s (top panel, solid line), and mass fraction of free neutrons (bottom 
panel), as functions of depth coordinate P/g. The composition from Haensel & Zdunik 1 2008 1 is shown (open circles) for comparison. 



Nuclear heating and neutrino cooling 

Following Brown ( 2000| l, we define a smooth heating distribution in the crust, rather than resolving the heating from individual 
reaction layers. We choose our heating function to be such that dL nuc /dln;y = const, and we do this separately in both the inner 



^nuc / 

crust, and in the outer crust where the pressure is P > 10 27 ergs cm 3 . The integrated nuclear luminosity is plotted in Fig. 16[ We 



normalized the heat distribution so that the total heat deposited, per accreted nucleon, into the inner crust is 1 .5 Me V (cf. Haensel 
& Zdunik|1990 2003 , 2008) and the total heat deposited, per accreted nucleon, into the outer crust is 0.2 MeV (cf. Gupta et al. 
2007) . 



For the neutrino cooling, our model includes (for a review of neutrino emission mechanisms, see Yakovlev et al. 2001 ) neutrino 
cooling from electron-nucleus bremsstrahlung. The neutrino emissivity from neutrons paired in th e state in the inner c rust is 
suppressed by a factor (vp/c) 2 (Leinson & Perez 2006[[Sedrakian et al.|2007 1. Recent calculations ( Steiner & Reddy|2009| show 
that this su ppression follows from conservation of baryon vector current. The pair, photo, and plasmon emissivities ( |Schinder| 
et al.|l9 87) do not contribute substantially at the temperatures of interest. 



Thermal conductivities 



Our implementation of the thermal conductivities mediated by electron-ion scattering follows that of Potekhin et al. ( 1999} and 
Gnedin et al. ( 2001| l. We compute the electron thermal transport in the relaxation-time approximation, 



K 



(A3) 



where m* = (p F lc + m e ) ' , with p F being the Fermi momentum, and v is the scattering frequency. In the ocean, v is set by 
electron-ion scattering. As the ions crystallize, electron-phonon scattering mediates the thermal transport. Where the temperature 
is above the Debye temperature, the scattering frequency is approximately 

k B T 
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Fig. 16. — Integrated nuclear heating, divided by the proper mass accretion rate, in the crust as a function of column. 



where a = e 2 /(hc) is the fine-structure constant. In the inner crust, the electron-ion scattering frequency is strongly reduced for 
T < T p , the plasma temperature, and impurity scattering becomes dominant with scattering frequency 



VeQ 



P 2 F v F 
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(A5) 



where pp and Vp are the momentum and velocity of electrons at the Fermi surface and the impurity parameter Q lmp = 
n ZL Hi n i(Zj - (Z)) 2 measures the distribution of nuclide charge numbers. Although we do include electron-electron scattering 



( Urpin & Yakovlev| 1980| |Potekhin et al.|1997 1 in our calculatio n, it does not affect the overall thermal conductivity. 

For the Coulomb logarithm term Ai mp we use the formula of Potekhin et al. ( 1999 ) with the modification that the structure factor 
is set to unity, reflecting the lack of long-range correlations in the impurities. With this modification Aj mp becomes (Potekhin, 
private communication) 
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where kp = pp/fr is the Fermi wavevector, /3 = v F /c is the Fermi velocity, and q s * k-pp is the Thomas-Fermi wavevector. As 
noted by |Brown et al.| | |2002") , this gives a result for the Coulomb logarithm that is similar to that proposed by |Itoh & Koh yama 
( |1993[ ). The two formulae agree if one makes the substitution kpa — > fc F /(fc XF /2), where a is the mean inter-nuclei spacing, in 
the fit by Itoh & Kohyama ( 1993 I. The Thomas-Fermi screening length exceeds the inter-nuclei spacing and produces a larger 
Ai m p (by about a factor of 2) and hence a lower thermal conductivity than does the fit by Itoh & Kohyama ( |1993 1. Molecular 
dynamics simulations ( Horowitz et al.|2008a l using a mixture of rp-process ashes ( Uupta et al.|2W/) , find that the impurities 
are not distributed uniformly. As a result, Horowitz et al. ( 2008a|l compute a thermal conductivity that is 30% lower than that 



computed assuming a one-component plasma and impurities computed with the static structure factor of Itoh & Kohyama ( 1993 1. 

Finally, Agui lera et al.| ( [2008] > recently calculated the thermal conductivity in the inner crust due to superfluid heat conduction. 
Their Figure 2 shows that, for T = 10 8 K, superfluid heat conduction is comparable to heat conduction by electrons for p n[ j < 
p S 2p n d, where p n d is the neutron drip density, and less important at lower temperatures. Therefore, we do not expect that our 
results would change significantly if superfluid heat conduction was included, unlike the case of magnetars, for which the strong 
magnetic field suppresses electron transport perpendicular to the magnetic field. 



